options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)

execute mcmc sampling

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

see mcmc result and parameters

seeMCMC=function(mcmc,exc='',ch=T){ # not see parameters str1..., str2,... using regex as exc='[str1,str2,...]' 
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(4,1),mar=c(1,5,1,1))
      drw[,1,] |> plot(type='l',xlab='',ylab=pr)
      drw[,2,] |> plot(type='l',xlab='',ylab=pr)
      drw[,3,] |> plot(type='l',xlab='',ylab=pr)
      drw[,4,] |> plot(type='l',xlab='',ylab=pr)
      par(mar=c(3,5,3,3))
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}
fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)

mcmc=goMCMC(mdl,data,smp,wrm)

mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')


y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')


m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')
}

hierarchical model

class c=1~k 
b0c~N(b00,sb0)  
b1c~N(b10,sb1)  
yc~N(b0c+b1c*x,s)
n=50
k=9
x=runif(n,0,20)
c=sample(letters[1:k],n,replace=T)
b00=rnorm(k,10,5)
b0=b00[factor(c)]
b10=rnorm(k,2,1)
b1=b10[factor(c)]
y=rnorm(n,b0+b1*x,5)
qplot(x,y,shape=c,size=I(2))+
  scale_shape_manual(values=1:k)

qplot(x,y,col=c)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

estimate as no class

ex8-0.stan

y~N(b00+b10*x,s)

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
mdl=cmdstan_model('./ex8-0.stan') 
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
fn(mdl,data)
## Initial log joint probability = -4279.57 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       21      -147.802    0.00185173    0.00427274      0.9623      0.9623       34    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.2 seconds.
##  variable estimate
##    lp__    -147.80
##    b0         8.30
##    b1         3.06
##    s         11.66
##    m0[1]     12.68
##    m0[2]     54.80
##    m0[3]     50.94
##    m0[4]     51.73
##    m0[5]     40.74
##    m0[6]     58.71
##    m0[7]     14.48
##    m0[8]     21.19
##    m0[9]     42.97
##    m0[10]    57.17
##    m0[11]    27.00
##    m0[12]    61.03
##    m0[13]    60.31
##    m0[14]    32.99
##    m0[15]    19.88
##    m0[16]    51.52
##    m0[17]    26.94
##    m0[18]    23.45
##    m0[19]    67.97
##    m0[20]    57.35
##    m0[21]    15.96
##    m0[22]    32.89
##    m0[23]    18.16
##    m0[24]    65.77
##    m0[25]    16.15
##    m0[26]    37.52
##    m0[27]    14.73
##    m0[28]    65.51
##    m0[29]    47.44
##    m0[30]    13.93
##    m0[31]    25.28
##    m0[32]    36.12
##    m0[33]    47.17
##    m0[34]    60.20
##    m0[35]    17.08
##    m0[36]    36.67
##    m0[37]    39.22
##    m0[38]    28.08
##    m0[39]     8.73
##    m0[40]    24.25
##    m0[41]    38.03
##    m0[42]    32.88
##    m0[43]    69.21
##    m0[44]    20.78
##    m0[45]    18.62
##    m0[46]    52.52
##    m0[47]    64.15
##    m0[48]    33.00
##    m0[49]    54.09
##    m0[50]    39.93
##    y0[1]     14.25
##    y0[2]     61.30
##    y0[3]     54.72
##    y0[4]     61.84
##    y0[5]     48.18
##    y0[6]     67.23
##    y0[7]      8.33
##    y0[8]     12.94
##    y0[9]     35.61
##    y0[10]    40.00
##    y0[11]    22.85
##    y0[12]    34.27
##    y0[13]    65.45
##    y0[14]    20.48
##    y0[15]    43.05
##    y0[16]    51.08
##    y0[17]    27.82
##    y0[18]    14.98
##    y0[19]    58.33
##    y0[20]    59.86
##    y0[21]     5.64
##    y0[22]    34.98
##    y0[23]    14.57
##    y0[24]    64.19
##    y0[25]     3.89
##    y0[26]    28.39
##    y0[27]    29.21
##    y0[28]    68.45
##    y0[29]    40.53
##    y0[30]    16.80
##    y0[31]    27.90
##    y0[32]    46.10
##    y0[33]    41.44
##    y0[34]    62.64
##    y0[35]    23.20
##    y0[36]    32.50
##    y0[37]    36.19
##    y0[38]    42.13
##    y0[39]     7.30
##    y0[40]    11.17
##    y0[41]    43.63
##    y0[42]    30.78
##    y0[43]    64.93
##    y0[44]    30.09
##    y0[45]    20.63
##    y0[46]    87.02
##    y0[47]    66.83
##    y0[48]    52.52
##    y0[49]    21.48
##    y0[50]    41.61
##    m1[1]      8.73
##    m1[2]     15.45
##    m1[3]     22.17
##    m1[4]     28.89
##    m1[5]     35.61
##    m1[6]     42.33
##    m1[7]     49.05
##    m1[8]     55.77
##    m1[9]     62.49
##    m1[10]    69.21
##    y1[1]    -18.31
##    y1[2]     31.64
##    y1[3]     29.22
##    y1[4]     19.51
##    y1[5]     44.08
##    y1[6]     42.64
##    y1[7]     55.84
##    y1[8]     53.35
##    y1[9]     59.47
##    y1[10]    74.10
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.4 seconds.
## 
##  variable    mean  median    sd   mad      q5     q95 rhat ess_bulk ess_tail
##    lp__   -146.88 -146.57  1.24  1.00 -149.47 -145.54 1.00      631      755
##    b0        8.45    8.50  3.46  3.53    2.80   13.90 1.01      279      265
##    b1        3.04    3.04  0.30  0.31    2.56    3.54 1.01      355      507
##    s        12.21   12.08  1.28  1.22   10.34   14.44 1.00     1923     1288
##    m0[1]    12.80   12.83  3.09  3.17    7.71   17.71 1.01      280      274
##    m0[2]    54.73   54.79  2.36  2.33   50.81   58.51 1.00     1351      974
##    m0[3]    50.89   50.93  2.11  2.05   47.34   54.25 1.00     1712     1126
##    m0[4]    51.68   51.73  2.16  2.09   48.02   55.16 1.00     1699     1086
##    m0[5]    40.74   40.71  1.73  1.68   37.82   43.49 1.00      930     1163
##    m0[6]    58.62   58.62  2.64  2.65   54.23   62.80 1.00     1068      928
##    m0[7]    14.60   14.64  2.94  3.01    9.82   19.29 1.01      281      288
##    m0[8]    21.28   21.25  2.43  2.47   17.33   25.16 1.01      293      355
##    m0[9]    42.96   42.97  1.77  1.71   40.02   45.79 1.00     1217     1198
##    m0[10]   57.09   57.12  2.52  2.52   52.92   61.10 1.00     1165      986
##    m0[11]   27.06   27.05  2.06  2.09   23.65   30.40 1.01      327      418
##    m0[12]   60.93   60.92  2.82  2.82   56.22   65.44 1.00      948      938
##    m0[13]   60.21   60.20  2.76  2.75   55.59   64.62 1.00      981      928
##    m0[14]   33.02   33.01  1.80  1.77   30.02   35.93 1.01      424      641
##    m0[15]   19.98   19.98  2.52  2.57   15.86   23.99 1.01      289      349
##    m0[16]   51.47   51.52  2.15  2.08   47.84   54.92 1.00     1703     1099
##    m0[17]   27.00   27.00  2.06  2.09   23.59   30.34 1.01      326      418
##    m0[18]   23.53   23.51  2.27  2.29   19.79   27.17 1.01      302      408
##    m0[19]   67.85   67.87  3.39  3.37   62.14   73.37 1.00      739      876
##    m0[20]   57.27   57.31  2.54  2.52   53.10   61.29 1.00     1153      986
##    m0[21]   16.07   16.12  2.82  2.92   11.48   20.58 1.01      282      313
##    m0[22]   32.92   32.91  1.80  1.78   29.91   35.83 1.01      421      626
##    m0[23]   18.26   18.28  2.65  2.72   13.95   22.46 1.01      285      334
##    m0[24]   65.65   65.69  3.21  3.20   60.32   70.78 1.00      791      874
##    m0[25]   16.26   16.30  2.81  2.90   11.69   20.75 1.01      283      313
##    m0[26]   37.53   37.51  1.71  1.67   34.66   40.30 1.01      638      970
##    m0[27]   14.85   14.89  2.92  3.00   10.11   19.52 1.01      282      288
##    m0[28]   65.40   65.44  3.18  3.17   60.11   70.48 1.00      797      874
##    m0[29]   47.41   47.45  1.93  1.90   44.10   50.49 1.00     1720     1265
##    m0[30]   14.05   14.09  2.98  3.05    9.17   18.80 1.01      281      281
##    m0[31]   25.34   25.34  2.16  2.18   21.79   28.83 1.01      312      421
##    m0[32]   36.14   36.10  1.73  1.70   33.24   38.92 1.01      545      930
##    m0[33]   47.14   47.17  1.92  1.89   43.86   50.19 1.00     1715     1287
##    m0[34]   60.10   60.09  2.75  2.75   55.49   64.50 1.00      985      928
##    m0[35]   17.18   17.22  2.73  2.80   12.69   21.52 1.01      284      311
##    m0[36]   36.68   36.64  1.72  1.69   33.76   39.46 1.01      578      940
##    m0[37]   39.22   39.20  1.71  1.67   36.32   41.97 1.00      782     1027
##    m0[38]   28.14   28.16  2.00  2.04   24.77   31.40 1.01      336      482
##    m0[39]    8.88    8.94  3.42  3.49    3.29   14.27 1.01      279      266
##    m0[40]   24.32   24.32  2.22  2.25   20.68   27.89 1.01      306      412
##    m0[41]   38.04   38.01  1.71  1.69   35.13   40.80 1.00      680      994
##    m0[42]   32.92   32.91  1.80  1.78   29.91   35.82 1.01      421      626
##    m0[43]   69.08   69.10  3.50  3.46   63.20   74.77 1.00      716      866
##    m0[44]   20.87   20.84  2.46  2.50   16.85   24.78 1.01      291      366
##    m0[45]   18.72   18.73  2.62  2.67   14.47   22.88 1.01      286      336
##    m0[46]   52.46   52.53  2.21  2.16   48.75   56.02 1.00     1678     1016
##    m0[47]   64.03   64.08  3.07  3.07   58.95   68.97 1.00      838      908
##    m0[48]   33.03   33.02  1.80  1.77   30.03   35.94 1.01      424      641
##    m0[49]   54.02   54.10  2.31  2.26   50.16   57.73 1.00     1459      947
##    m0[50]   39.93   39.92  1.72  1.67   37.03   42.68 1.00      847     1141
##    y0[1]    13.35   13.54 12.52 12.37   -7.57   33.51 1.00     1449     1805
##    y0[2]    55.12   55.08 12.38 12.19   35.24   75.10 1.00     2165     2105
##    y0[3]    50.92   50.63 12.63 12.04   30.66   72.00 1.00     1995     1856
##    y0[4]    51.53   51.47 12.65 12.50   31.08   72.46 1.00     1869     1804
##    y0[5]    40.84   40.52 12.56 12.34   20.60   61.41 1.00     2014     1931
##    y0[6]    58.49   58.73 12.22 11.92   38.19   78.12 1.00     2011     1906
##    y0[7]    14.26   14.56 12.60 12.81   -6.51   35.05 1.00     1548     1931
##    y0[8]    21.31   21.15 12.53 11.96    0.18   41.57 1.00     1797     1972
##    y0[9]    43.31   43.28 12.43 12.26   23.21   63.93 1.00     1929     1798
##    y0[10]   56.95   57.38 12.22 12.20   36.23   76.93 1.00     1911     1973
##    y0[11]   26.92   26.68 12.37 12.01    6.66   47.09 1.00     2004     1890
##    y0[12]   60.91   60.98 12.25 11.90   40.97   80.92 1.00     1843     1900
##    y0[13]   60.49   60.43 12.49 12.34   40.02   81.11 1.00     1990     2058
##    y0[14]   33.42   33.81 12.22 12.14   12.51   53.44 1.00     1959     1996
##    y0[15]   20.25   20.89 12.50 12.76   -0.67   40.87 1.00     1998     1828
##    y0[16]   51.24   51.16 12.41 12.30   31.41   71.70 1.00     1650     2008
##    y0[17]   27.07   27.03 12.21 12.44    6.68   46.84 1.00     2015     1933
##    y0[18]   23.47   23.62 12.50 12.26    2.88   43.42 1.00     2011     1934
##    y0[19]   67.53   67.60 12.95 12.95   46.69   88.10 1.00     1919     1972
##    y0[20]   57.14   57.31 12.80 12.66   36.19   77.47 1.00     1773     1969
##    y0[21]   15.67   16.04 12.69 12.35   -5.39   36.31 1.00     1587     1757
##    y0[22]   33.13   32.94 12.31 12.14   12.82   53.35 1.00     2131     2082
##    y0[23]   18.27   18.24 12.13 11.75   -1.72   38.37 1.00     1910     1779
##    y0[24]   66.09   65.73 12.83 12.57   44.54   87.55 1.00     1810     1795
##    y0[25]   15.95   15.94 12.26 11.73   -4.13   36.33 1.00     1679     1918
##    y0[26]   37.10   36.98 12.74 12.05   15.71   58.49 1.00     1753     1882
##    y0[27]   15.13   14.93 12.64 12.41   -5.68   35.15 1.00     1733     1741
##    y0[28]   65.16   65.09 12.80 12.13   43.79   86.44 1.00     1744     1910
##    y0[29]   47.08   46.91 12.97 12.94   25.67   68.26 1.00     2042     1779
##    y0[30]   13.75   14.07 12.73 12.92   -7.48   34.61 1.00     1816     1766
##    y0[31]   25.29   25.18 12.70 12.51    4.56   46.49 1.00     1719     1709
##    y0[32]   36.55   36.79 12.61 12.82   15.45   57.82 1.00     1927     2057
##    y0[33]   47.02   47.34 12.52 12.92   26.37   67.02 1.00     2116     1882
##    y0[34]   60.43   60.45 12.58 12.39   39.31   80.96 1.00     1845     1907
##    y0[35]   17.17   16.99 12.48 12.15   -3.18   37.83 1.00     1899     1929
##    y0[36]   36.49   36.58 12.06 12.12   16.42   56.75 1.00     1936     1659
##    y0[37]   39.27   39.16 12.21 11.99   19.05   58.70 1.00     2013     1818
##    y0[38]   28.09   28.56 12.18 12.22    7.97   47.40 1.00     1911     1923
##    y0[39]    9.07    8.89 12.96 12.50  -11.68   30.28 1.00     1314     1725
##    y0[40]   23.97   23.87 12.62 12.00    3.71   45.07 1.00     1743     1719
##    y0[41]   38.17   38.15 12.19 12.25   18.53   58.39 1.00     1810     1812
##    y0[42]   33.04   32.97 12.21 12.35   13.06   53.00 1.00     1815     1794
##    y0[43]   68.92   68.78 13.23 13.36   47.46   90.55 1.00     1976     1884
##    y0[44]   20.60   20.62 12.49 12.24    0.27   41.04 1.00     1819     1932
##    y0[45]   18.07   18.07 12.56 12.48   -2.96   38.29 1.00     1934     1861
##    y0[46]   52.65   52.69 12.26 12.83   33.12   72.39 1.00     1888     1971
##    y0[47]   64.00   63.88 12.82 12.63   43.21   85.26 1.00     1972     1946
##    y0[48]   32.88   32.61 12.08 11.71   12.92   52.86 1.00     2129     1933
##    y0[49]   53.50   53.88 12.39 12.35   32.66   72.97 1.00     1986     2014
##    y0[50]   40.28   40.31 12.34 12.43   20.42   60.31 1.00     1705     1838
##    m1[1]     8.88    8.94  3.42  3.49    3.29   14.27 1.01      279      266
##    m1[2]    15.57   15.63  2.86  2.95   10.92   20.15 1.01      282      307
##    m1[3]    22.26   22.23  2.36  2.39   18.40   26.03 1.01      296      369
##    m1[4]    28.94   28.96  1.96  1.99   25.65   32.16 1.01      345      528
##    m1[5]    35.63   35.59  1.74  1.70   32.72   38.44 1.01      517      857
##    m1[6]    42.32   42.32  1.75  1.70   39.39   45.15 1.00     1126     1224
##    m1[7]    49.01   49.05  2.01  1.96   45.58   52.20 1.00     1725     1196
##    m1[8]    55.70   55.74  2.42  2.40   51.67   59.58 1.00     1264      974
##    m1[9]    62.39   62.42  2.93  2.93   57.49   67.11 1.00      890      920
##    m1[10]   69.08   69.10  3.50  3.46   63.20   74.77 1.00      716      866
##    y1[1]     8.82    9.13 12.72 12.63  -11.69   29.96 1.00     1344     1656
##    y1[2]    15.81   15.99 12.51 11.84   -5.60   36.09 1.00     1735     1715
##    y1[3]    22.23   21.97 12.49 12.22    1.96   43.32 1.00     1841     1683
##    y1[4]    28.70   28.78 12.33 12.30    8.18   48.76 1.00     1565     2015
##    y1[5]    36.00   35.29 12.12 11.85   17.15   56.29 1.00     1625     1792
##    y1[6]    42.29   42.36 12.37 12.33   22.01   61.89 1.00     1852     1927
##    y1[7]    49.33   49.74 12.21 12.26   28.18   68.34 1.00     1892     1969
##    y1[8]    55.34   54.93 12.65 12.52   34.41   76.11 1.00     1942     1740
##    y1[9]    62.55   62.50 12.80 12.68   41.69   83.36 1.00     2011     1860
##    y1[10]   68.85   69.12 12.68 12.62   48.24   89.15 1.00     1834     1786

estimate as independent class

all b0l,b1l do not have a distribution  
class c=1~k 

ex12-1.stan

yc~N(b0c+b1c*x,s)

data{
  int N;
  int K;
  vector[N] x;
  vector[N] y;
  array[N] int c;
}
parameters{
  vector[K] b0;
  vector[K] b1;
  real<lower=0> s;
}
model{
  for(i in 1:N)
    y[i]~normal(b0[c[i]]+b1[c[i]]*x[i],s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0[c[i]]+b1[c[i]]*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
}
mdl=cmdstan_model('./ex12-1.stan') 
data=list(N=n,K=k,x=x,y=y,c=factor(c))

mle=mdl$optimize(data=data)
## Initial log joint probability = -587644 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmpwb3xsj/model-20a93d7d50fa.stan', line 15, column 4 to column 42) 
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmpwb3xsj/model-20a93d7d50fa.stan', line 15, column 4 to column 42) 
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmpwb3xsj/model-20a93d7d50fa.stan', line 15, column 4 to column 42) 
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
##       99      -86.3037      0.232138      0.816477           1           1      154    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      199      -85.7831     0.0876634      0.790418      0.2119           1      265    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      299      -85.4373    0.00807286     0.0674109      0.9062      0.9062      370    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      399      -85.4126    0.00123051    0.00769343           1           1      481    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      499      -85.4118   0.000725348    0.00279152           1           1      587    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -85.41
##    b0[1]     -2.33
##    b0[2]     -0.42
##    b0[3]     26.04
##    b0[4]      5.46
##    b0[5]     16.22
##    b0[6]      1.71
##    b0[7]     12.72
##    b0[8]     11.71
##    b0[9]     13.13
##    b1[1]      0.87
##    b1[2]      4.00
##    b1[3]      3.55
##    b1[4]      3.99
##    b1[5]      3.19
##    b1[6]      3.96
##    b1[7]      2.55
##    b1[8]      2.21
##    b1[9]      2.89
##    s          3.35
##    m0[1]     17.27
##    m0[2]     51.53
##    m0[3]     53.47
##    m0[4]     43.07
##    m0[5]      6.88
##    m0[6]     11.98
##    m0[7]     22.66
##    m0[8]     16.45
##    m0[9]     36.75
##    m0[10]    67.15
##    m0[11]    24.04
##    m0[12]    74.32
##    m0[13]    70.42
##    m0[14]    33.66
##    m0[15]    24.09
##    m0[16]    76.26
##    m0[17]    29.80
##    m0[18]    25.36
##    m0[19]    69.58
##    m0[20]    53.65
##    m0[21]    11.62
##    m0[22]    36.39
##    m0[23]    22.46
##    m0[24]    67.50
##    m0[25]    15.71
##    m0[26]    37.81
##    m0[27]    18.08
##    m0[28]    53.02
##    m0[29]    39.97
##    m0[30]    -0.73
##    m0[31]    29.19
##    m0[32]    35.98
##    m0[33]    50.44
##    m0[34]    49.18
##    m0[35]    20.04
##    m0[36]    39.97
##    m0[37]    40.03
##    m0[38]    27.31
##    m0[39]    13.54
##    m0[40]    26.03
##    m0[41]    37.53
##    m0[42]    29.46
##    m0[43]    70.75
##    m0[44]    17.86
##    m0[45]    18.94
##    m0[46]    62.30
##    m0[47]    90.92
##    m0[48]    33.33
##    m0[49]    63.94
##    m0[50]    40.97
##    y0[1]     15.37
##    y0[2]     44.33
##    y0[3]     54.84
##    y0[4]     46.67
##    y0[5]     -0.74
##    y0[6]     14.73
##    y0[7]     25.49
##    y0[8]     14.23
##    y0[9]     28.75
##    y0[10]    72.09
##    y0[11]    21.59
##    y0[12]    76.05
##    y0[13]    69.59
##    y0[14]    31.42
##    y0[15]    28.80
##    y0[16]    77.22
##    y0[17]    29.13
##    y0[18]    31.56
##    y0[19]    67.90
##    y0[20]    56.16
##    y0[21]    10.52
##    y0[22]    34.49
##    y0[23]    22.78
##    y0[24]    70.25
##    y0[25]    13.85
##    y0[26]    37.22
##    y0[27]    22.34
##    y0[28]    54.94
##    y0[29]    41.55
##    y0[30]     4.06
##    y0[31]    28.63
##    y0[32]    35.30
##    y0[33]    45.37
##    y0[34]    45.88
##    y0[35]    24.26
##    y0[36]    40.22
##    y0[37]    42.55
##    y0[38]    24.25
##    y0[39]    15.27
##    y0[40]    26.94
##    y0[41]    36.98
##    y0[42]    24.45
##    y0[43]    70.04
##    y0[44]    16.32
##    y0[45]    12.37
##    y0[46]    60.08
##    y0[47]    89.04
##    y0[48]    36.06
##    y0[49]    60.55
##    y0[50]    34.71
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 finished in 1.7 seconds.
## Chain 3 finished in 1.8 seconds.
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 1.9 seconds.
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 finished in 2.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 1.8 seconds.
## Total execution time: 2.1 seconds.
seeMCMC(mcmc,'m')
##  variable   mean median    sd   mad      q5    q95 rhat ess_bulk ess_tail
##    lp__   -96.02 -95.58  3.74  3.79 -102.85 -90.47 1.01      580      988
##    b0[1]   -2.24  -2.18  4.61  4.52   -9.67   5.26 1.00     2567     1463
##    b0[2]   -0.22  -0.21  5.99  5.84  -10.06   9.45 1.01     1968     1154
##    b0[3]   27.80  26.24 24.69 24.65  -10.91  70.14 1.01      271      363
##    b0[4]    5.51   5.55  3.38  3.20   -0.03  10.99 1.01     3016     1517
##    b0[5]   16.23  16.19  5.04  5.12    7.81  24.51 1.01     2161     1363
##    b0[6]    1.82   1.83  5.59  5.61   -7.43  11.23 1.00     1842     1302
##    b0[7]   12.71  12.71  2.94  2.87    7.82  17.55 1.00     2421     1149
##    b0[8]   11.43  11.63  7.24  7.21   -0.76  22.69 1.00     1331     1016
##    b0[9]   13.18  13.23  2.08  2.09    9.79  16.61 1.00     2389     1199
##    b1[1]    0.85   0.85  0.40  0.40    0.18   1.49 1.00     2459     1664
##    b1[2]    3.99   3.98  0.64  0.62    2.96   5.06 1.01     1986     1096
##    b1[3]    3.44   3.54  1.52  1.56    0.84   5.83 1.01      275      416
##    b1[4]    3.99   3.99  0.37  0.36    3.40   4.60 1.00     2744     1419
##    b1[5]    3.19   3.19  0.36  0.36    2.60   3.80 1.01     2145     1490
##    b1[6]    3.95   3.93  0.97  0.94    2.36   5.56 1.00     1804     1388
##    b1[7]    2.56   2.56  0.31  0.30    2.04   3.07 1.00     2506     1344
##    b1[8]    2.23   2.22  0.51  0.52    1.44   3.10 1.00     1360     1081
##    b1[9]    2.89   2.88  0.18  0.18    2.60   3.18 1.00     2231     1580
##    s        4.33   4.27  0.56  0.55    3.53   5.35 1.00     1182     1456
##    m0[1]   17.31  17.33  1.89  1.88   14.22  20.43 1.00     2370     1203
##    m0[2]   51.61  51.64  2.73  2.59   46.99  56.14 1.00     2189     1847
##    m0[3]   53.46  53.45  1.53  1.49   51.00  55.91 1.00     2175     1338
##    m0[4]   43.15  43.18  1.80  1.75   40.09  46.04 1.01     2108     1288
##    m0[5]    6.78   6.75  2.58  2.51    2.49  10.89 1.00     1923     1557
##    m0[6]   11.78  11.84  3.79  3.78    5.72  17.81 1.00     2001     1628
##    m0[7]   22.68  22.64  4.38  4.41   15.31  29.80 1.01     2151     1382
##    m0[8]   16.61  16.58  3.46  3.34   10.98  22.20 1.00     1945     1250
##    m0[9]   36.75  36.83  2.16  2.06   33.06  40.23 1.00     1590     1646
##    m0[10]  67.22  67.21  2.18  2.14   63.68  70.87 1.01     1862     1209
##    m0[11]  24.18  24.13  2.46  2.40   20.16  28.10 1.00     1956     1358
##    m0[12]  74.38  74.33  4.34  4.39   67.44  81.31 1.00     2219     1529
##    m0[13]  70.49  70.45  2.37  2.36   66.63  74.41 1.01     1875     1366
##    m0[14]  33.71  33.59  3.43  3.40   28.10  39.34 1.00     1799     1249
##    m0[15]  24.12  24.16  1.60  1.59   21.42  26.72 1.01     2323     1377
##    m0[16]  76.42  76.45  4.29  4.16   69.14  83.46 1.00      402      654
##    m0[17]  29.85  29.86  2.19  2.11   26.27  33.53 1.00     2267     1784
##    m0[18]  25.38  25.37  1.81  1.73   22.45  28.48 1.00     2515     1494
##    m0[19]  69.54  69.50  2.24  2.17   65.97  73.26 1.00     2134     1482
##    m0[20]  53.75  53.76  2.94  2.80   48.77  58.64 1.00     2228     1831
##    m0[21]  11.71  11.76  3.48  3.48    5.98  17.50 1.00     1875     1438
##    m0[22]  36.40  36.39  1.29  1.31   34.25  38.38 1.00     2193     1633
##    m0[23]  22.49  22.53  1.66  1.67   19.70  25.20 1.00     2339     1320
##    m0[24]  67.46  67.43  2.13  2.06   64.04  70.99 1.00     2240     1481
##    m0[25]  15.77  15.74  2.72  2.57   11.36  20.36 1.00     2849     1555
##    m0[26]  37.92  37.89  1.69  1.64   35.13  40.64 1.00     2158     1661
##    m0[27]  18.09  18.06  2.40  2.31   14.17  22.06 1.00     2443     1154
##    m0[28]  53.21  53.22  3.11  3.10   48.30  58.21 1.00     1729     1347
##    m0[29]  40.01  40.04  1.84  1.77   36.88  43.00 1.00     1872     1657
##    m0[30]  -0.67  -0.61  4.01  3.95   -7.09   5.98 1.00     2514     1424
##    m0[31]  29.21  29.26  1.43  1.45   26.82  31.50 1.01     2268     1649
##    m0[32]  36.09  36.09  1.64  1.62   33.38  38.76 1.00     2129     1624
##    m0[33]  50.51  50.43  2.92  2.97   45.79  55.21 1.00     2108     1593
##    m0[34]  49.33  49.31  2.43  2.43   45.41  53.31 1.00     1888     1548
##    m0[35]  20.05  20.07  2.22  2.16   16.41  23.78 1.00     2433     1171
##    m0[36]  39.98  39.97  1.28  1.27   37.91  42.00 1.00     2162     1671
##    m0[37]  40.13  40.09  1.80  1.74   37.14  43.06 1.00     2181     1559
##    m0[38]  27.37  27.38  2.42  2.44   23.33  31.29 1.00     1830     1074
##    m0[39]  13.59  13.64  2.06  2.07   10.22  16.99 1.00     2385     1199
##    m0[40]  26.05  26.04  1.77  1.70   23.20  29.07 1.00     2493     1479
##    m0[41]  37.58  37.56  1.63  1.58   34.89  40.31 1.00     2071     1571
##    m0[42]  29.38  29.53  3.41  3.37   23.50  34.73 1.00     1322     1022
##    m0[43]  70.71  70.67  2.30  2.24   67.02  74.48 1.00     2132     1547
##    m0[44]  17.94  17.91  2.46  2.49   13.96  21.99 1.00     1915     1619
##    m0[45]  19.00  18.97  2.55  2.43   14.93  23.26 1.00     2751     1583
##    m0[46]  62.36  62.39  1.98  1.98   59.10  65.55 1.01     1892     1194
##    m0[47]  90.62  90.61  4.53  4.42   83.09  97.87 1.01      570      816
##    m0[48]  33.37  33.32  1.54  1.47   30.83  36.06 1.00     2226     1576
##    m0[49]  64.00  64.04  2.03  2.02   60.64  67.35 1.01     1865     1194
##    m0[50]  41.07  41.00  1.87  1.83   38.04  44.12 1.00     2180     1557
##    y0[1]   17.19  17.12  4.74  4.59    9.51  24.99 1.00     1836     1705
##    y0[2]   51.59  51.47  5.08  5.09   43.38  59.91 1.00     2035     1890
##    y0[3]   53.53  53.66  4.53  4.50   46.42  60.91 1.00     2035     1972
##    y0[4]   43.05  43.16  4.62  4.59   35.26  50.20 1.00     1794     1974
##    y0[5]    6.65   6.52  5.18  4.84   -1.85  15.45 1.00     1831     1714
##    y0[6]   11.86  12.01  5.65  5.40    2.80  20.95 1.00     1873     1374
##    y0[7]   22.67  22.64  6.05  5.87   12.88  32.87 1.00     2248     1767
##    y0[8]   16.58  16.69  5.59  5.57    7.45  25.68 1.00     2001     1963
##    y0[9]   36.90  36.90  4.79  4.83   29.16  44.62 1.00     1779     1842
##    y0[10]  67.20  67.14  4.91  4.97   59.21  75.29 1.00     1902     1951
##    y0[11]  24.05  24.01  4.86  4.77   16.20  31.84 1.00     2004     1752
##    y0[12]  74.36  74.32  6.13  6.10   64.45  84.38 1.00     1931     1931
##    y0[13]  70.35  70.23  5.01  4.85   62.31  78.61 1.00     1991     1930
##    y0[14]  33.85  33.98  5.54  5.61   24.82  42.89 1.00     1973     1957
##    y0[15]  24.35  24.37  4.56  4.54   16.80  31.99 1.00     2009     1931
##    y0[16]  76.46  76.41  6.19  6.09   66.71  87.00 1.00      760     1472
##    y0[17]  29.82  29.82  4.80  4.65   21.88  37.81 1.00     1962     1931
##    y0[18]  25.20  25.30  4.58  4.48   17.42  32.69 1.00     2036     1888
##    y0[19]  69.53  69.64  4.80  4.87   61.48  76.99 1.00     1978     1760
##    y0[20]  53.72  53.69  5.03  5.06   45.46  62.00 1.00     1999     1786
##    y0[21]  11.62  11.57  5.59  5.12    2.45  20.92 1.00     1997     1796
##    y0[22]  36.43  36.42  4.61  4.56   28.97  43.95 1.00     2083     2007
##    y0[23]  22.59  22.54  4.68  4.67   15.05  30.25 1.00     2074     1861
##    y0[24]  67.25  67.33  4.88  4.65   59.21  75.30 1.00     2118     2105
##    y0[25]  15.76  15.77  4.97  4.95    7.79  23.86 1.00     2073     1881
##    y0[26]  37.80  37.79  4.65  4.64   30.04  45.62 1.00     2079     2051
##    y0[27]  18.12  18.06  4.97  4.94    9.73  26.24 1.00     2156     1873
##    y0[28]  53.22  53.05  5.25  5.11   44.86  61.98 1.00     1922     1979
##    y0[29]  40.10  40.08  4.61  4.57   32.28  47.39 1.00     1973     1943
##    y0[30]  -0.77  -0.76  5.91  5.89  -10.30   8.98 1.00     2112     1842
##    y0[31]  29.13  29.23  4.60  4.54   21.48  36.72 1.00     1779     1825
##    y0[32]  36.10  36.15  4.61  4.57   28.62  43.64 1.00     2013     1898
##    y0[33]  50.68  50.90  5.18  5.22   42.23  58.68 1.00     2032     1926
##    y0[34]  49.35  49.20  5.00  4.77   41.11  57.78 1.00     2221     1931
##    y0[35]  20.25  20.20  5.04  4.84   12.09  28.56 1.00     2146     2056
##    y0[36]  39.92  39.98  4.68  4.43   32.38  47.64 1.00     2083     1933
##    y0[37]  40.15  40.06  4.74  4.75   32.56  48.12 1.00     2138     1816
##    y0[38]  27.42  27.53  4.88  4.83   19.46  35.16 1.00     1559     1785
##    y0[39]  13.49  13.56  4.91  4.71    5.39  21.43 1.00     2039     1955
##    y0[40]  26.01  25.92  4.73  4.69   18.39  33.69 1.00     2104     1913
##    y0[41]  37.63  37.81  4.72  4.70   29.37  45.38 1.00     2003     1695
##    y0[42]  29.18  29.12  5.61  5.39   19.85  38.20 1.00     1623     1678
##    y0[43]  70.80  70.82  4.85  4.81   62.66  79.01 1.00     1997     1892
##    y0[44]  17.95  17.93  5.16  5.07    9.59  26.45 1.00     1997     1781
##    y0[45]  18.91  18.96  5.11  5.04   10.39  27.44 1.00     2093     1767
##    y0[46]  62.35  62.36  4.74  4.52   54.37  70.01 1.00     2011     1970
##    y0[47]  90.65  90.71  6.28  6.26   80.12 100.80 1.01      855     1486
##    y0[48]  33.27  33.32  4.67  4.65   25.54  40.77 1.00     1894     1881
##    y0[49]  64.04  63.93  4.74  4.62   56.32  72.11 1.00     2022     1810
##    y0[50]  40.83  40.91  4.91  4.93   32.80  48.77 1.00     1988     1846

m0=mcmc$draws('m0')
smm0=summary(m0)

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

xy=tibble(x=x,c=c,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)

qplot(x,y,shape=c,size=I(2),col=I('red'))+
  scale_shape_manual(values=1:k) + 
  #geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  #geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  #geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  #geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

qplot(x,y,facets=~c,size=I(2),col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

estimate as class have relation

all b0l,b1l have a distribution  
class c=1~k 
b0c~N(b00,sb0)  
b1c~N(b10,sb1)  
yc~N(b0c+b1c*x,s)

ex12-2.stan

class have relation

data{
  int N;
  int K;
  vector[N] x;
  vector[N] y;
  array[N] int c;
}
parameters{
  real b00;
  real<lower=0,upper=100> sb0;
  vector[K] b0;
  real b10;
  real<lower=0,upper=100> sb1;
  vector[K] b1;
  real<lower=0,upper=100> s;
}
model{
  b0~normal(b00,sb0);
  b1~normal(b10,sb1);
  for(i in 1:N)
    y[i]~normal(b0[c[i]]+b1[c[i]]*x[i],s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0[c[i]]+b1[c[i]]*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
}
mdl=cmdstan_model('./ex12-2.stan') 
data=list(N=n,K=k,x=x,y=y,c=factor(c))

mle=mdl$optimize(data=data)
## Initial log joint probability = -550.305 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       99      -58.4879   0.000681448        169.74           1           1      133    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      199       34.9669     0.0759356   2.42048e+12      0.4956      0.4956      308    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      237       88.3164     0.0728787   1.11027e+16       1e-12       0.001      434  LS failed, Hessian reset  
## Finished in  0.1 seconds.
try(print(mle))
## Error : Fitting failed. Unable to print.
mcmc=goMCMC(mdl,data,wrm=1000,smp=2000)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 2 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 3 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 4 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 1 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 1 Iteration: 1001 / 3000 [ 33%]  (Sampling) 
## Chain 2 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 2 Iteration: 1001 / 3000 [ 33%]  (Sampling) 
## Chain 3 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 3 Iteration: 1001 / 3000 [ 33%]  (Sampling) 
## Chain 4 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 4 Iteration: 1001 / 3000 [ 33%]  (Sampling) 
## Chain 2 Iteration: 2000 / 3000 [ 66%]  (Sampling) 
## Chain 1 Iteration: 2000 / 3000 [ 66%]  (Sampling) 
## Chain 3 Iteration: 2000 / 3000 [ 66%]  (Sampling) 
## Chain 4 Iteration: 2000 / 3000 [ 66%]  (Sampling) 
## Chain 2 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 2 finished in 1.3 seconds.
## Chain 1 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 3 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 4 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 1 finished in 1.5 seconds.
## Chain 3 finished in 1.4 seconds.
## Chain 4 finished in 1.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 1.4 seconds.
## Total execution time: 1.5 seconds.
seeMCMC(mcmc,'m')
##  variable    mean  median   sd  mad      q5     q95 rhat ess_bulk ess_tail
##    lp__   -120.13 -119.90 4.70 4.33 -127.97 -113.43 1.01      320       77
##    b00       8.47    8.42 3.21 2.76    3.37   13.55 1.00     6940     5384
##    sb0       7.12    6.54 3.78 3.07    2.17   14.02 1.02      286       69
##    b0[1]     0.05    0.01 4.94 5.05   -7.99    8.11 1.01      853     1565
##    b0[2]     4.44    4.61 4.43 4.45   -3.24   11.26 1.00     4331     5670
##    b0[3]    12.90   11.80 8.00 6.72    1.64   27.29 1.00     6395     4832
##    b0[4]     6.78    6.91 3.02 2.91    1.71   11.64 1.00     7355     6071
##    b0[5]    13.35   13.14 4.19 4.25    6.91   20.46 1.00     1210     5644
##    b0[6]     5.54    5.78 3.89 3.79   -1.21   11.70 1.00     6341     6273
##    b0[7]    11.64   11.64 2.59 2.60    7.48   15.90 1.00     1307     5860
##    b0[8]     9.12    9.02 4.74 4.46    1.38   17.07 1.00     7109     5692
##    b0[9]    12.38   12.43 2.18 2.15    8.72   15.95 1.00      784      883
##    b10       3.02    3.02 0.48 0.45    2.25    3.78 1.00     5362     5121
##    sb1       1.31    1.23 0.45 0.37    0.76    2.12 1.00     5395     4455
##    b1[1]     0.77    0.77 0.43 0.44    0.05    1.48 1.01      740      794
##    b1[2]     3.49    3.47 0.48 0.47    2.72    4.31 1.00     4562     5129
##    b1[3]     4.32    4.37 0.51 0.46    3.39    5.07 1.00     6072     4744
##    b1[4]     3.85    3.86 0.34 0.33    3.30    4.41 1.00     7084     6019
##    b1[5]     3.37    3.38 0.31 0.31    2.86    3.85 1.00     1303     3245
##    b1[6]     3.30    3.28 0.69 0.68    2.22    4.46 1.00     6042     5911
##    b1[7]     2.66    2.65 0.28 0.28    2.20    3.11 1.00     1540     3804
##    b1[8]     2.40    2.41 0.35 0.34    1.83    2.95 1.00     5427     5524
##    b1[9]     2.95    2.94 0.18 0.18    2.65    3.24 1.00      967     2826
##    s         4.30    4.24 0.55 0.53    3.49    5.29 1.00     5213     5021
##    m0[1]    16.60   16.64 1.98 1.95   13.32   19.85 1.00      818     1484
##    m0[2]    52.04   52.02 2.59 2.55   47.82   56.30 1.00     6158     4052
##    m0[3]    53.45   53.47 1.56 1.48   50.88   55.96 1.00     8672     5249
##    m0[4]    43.15   43.14 1.76 1.69   40.28   46.09 1.00     2092      638
##    m0[5]     8.20    8.14 2.55 2.45    4.06   12.47 1.00     7573     5451
##    m0[6]    12.72   12.65 3.85 3.87    6.47   19.13 1.00     2292     1526
##    m0[7]    20.17   19.99 3.65 3.74   14.62   26.44 1.00     1316     5540
##    m0[8]    19.15   19.23 2.63 2.63   14.66   23.29 1.00     4307     6104
##    m0[9]    36.28   36.28 1.75 1.72   33.39   39.14 1.00     4419     5294
##    m0[10]   67.24   67.24 2.18 2.13   63.72   70.84 1.00     7276     5694
##    m0[11]   25.77   25.81 1.97 1.97   22.49   28.93 1.00     4885     6390
##    m0[12]   73.23   73.29 4.20 4.11   66.21   80.05 1.00     3895     3005
##    m0[13]   70.71   70.68 2.36 2.33   66.89   74.58 1.00     7049     5901
##    m0[14]   32.20   32.18 3.03 2.96   27.28   37.10 1.00     7700     5950
##    m0[15]   23.54   23.55 1.68 1.65   20.79   26.25 1.00      959     1928
##    m0[16]   73.95   73.93 2.89 2.86   69.28   78.69 1.00     6861     5956
##    m0[17]   30.27   30.28 2.12 2.09   26.78   33.73 1.00     2790     4921
##    m0[18]   24.80   24.78 1.66 1.66   22.16   27.53 1.00     3415     6575
##    m0[19]   69.85   69.85 2.28 2.18   66.14   73.56 1.00     5340     6021
##    m0[20]   54.25   54.22 2.78 2.78   49.69   58.82 1.00     5758     4152
##    m0[21]   13.81   13.87 2.59 2.52    9.45   17.96 1.00     7321     6577
##    m0[22]   36.06   36.07 1.34 1.34   33.88   38.24 1.00     2492     4813
##    m0[23]   21.88   21.90 1.75 1.72   19.01   24.71 1.00      917     1471
##    m0[24]   67.73   67.74 2.17 2.07   64.19   71.27 1.00     5682     6109
##    m0[25]   16.67   16.74 2.48 2.39   12.48   20.70 1.00     6312     5615
##    m0[26]   37.77   37.76 1.65 1.60   35.06   40.46 1.00     7586     5255
##    m0[27]   17.22   17.21 2.14 2.14   13.80   20.74 1.00     1580     6154
##    m0[28]   53.94   53.92 2.68 2.61   49.58   58.41 1.00     2227      641
##    m0[29]   39.78   39.78 1.68 1.65   36.99   42.58 1.00     2706     1312
##    m0[30]    1.46    1.42 4.27 4.34   -5.49    8.37 1.01      960     1767
##    m0[31]   28.73   28.73 1.50 1.49   26.31   31.16 1.00     1179     4341
##    m0[32]   36.17   36.18 1.60 1.56   33.56   38.79 1.00     7522     5110
##    m0[33]   48.78   48.77 2.56 2.49   44.65   52.99 1.00     6446     4969
##    m0[34]   49.78   49.75 2.24 2.16   46.13   53.51 1.00     2006      633
##    m0[35]   19.26   19.25 1.99 2.00   16.07   22.54 1.00     1821     6554
##    m0[36]   39.70   39.69 1.32 1.30   37.54   41.82 1.00     4850     4724
##    m0[37]   39.70   39.69 1.75 1.70   36.82   42.57 1.00     7524     5419
##    m0[38]   26.90   26.88 2.34 2.29   23.11   30.75 1.00     8121     5468
##    m0[39]   12.80   12.85 2.16 2.14    9.18   16.33 1.00      784      866
##    m0[40]   25.49   25.47 1.63 1.61   22.91   28.17 1.00     3954     6891
##    m0[41]   37.47   37.46 1.60 1.52   34.86   40.11 1.00     9189     5003
##    m0[42]   28.38   28.35 2.35 2.26   24.58   32.18 1.00     7522     6148
##    m0[43]   71.05   71.04 2.34 2.25   67.23   74.86 1.00     5171     6182
##    m0[44]   19.01   19.08 2.10 2.05   15.54   22.50 1.00     7791     5660
##    m0[45]   19.79   19.84 2.35 2.29   15.84   23.62 1.00     5408     6091
##    m0[46]   62.11   62.12 1.99 1.94   58.92   65.39 1.00     7457     5149
##    m0[47]   91.78   91.79 3.47 3.46   85.99   97.46 1.00     6769     6210
##    m0[48]   33.10   33.07 1.49 1.42   30.66   35.56 1.00     9143     5179
##    m0[49]   63.85   63.86 2.04 2.00   60.57   67.22 1.00     7425     5261
##    m0[50]   40.52   40.51 1.80 1.75   37.55   43.49 1.00     7474     5068
##    y0[1]    16.71   16.74 4.75 4.62    8.77   24.52 1.00     6739     7967
##    y0[2]    52.07   52.04 5.02 4.85   43.87   60.32 1.00     7624     7790
##    y0[3]    53.43   53.44 4.61 4.53   45.90   60.91 1.00     8458     7669
##    y0[4]    43.15   43.17 4.60 4.62   35.63   50.70 1.00     6986     7610
##    y0[5]     8.17    8.12 5.02 4.99    0.01   16.46 1.00     7784     7191
##    y0[6]    12.69   12.70 5.76 5.68    3.22   22.02 1.00     6821     7247
##    y0[7]    20.17   20.11 5.67 5.61   10.93   29.43 1.00     4865     6320
##    y0[8]    19.19   19.21 5.07 4.96   10.70   27.59 1.00     6405     6878
##    y0[9]    36.29   36.32 4.67 4.47   28.51   43.90 1.00     8123     7924
##    y0[10]   67.26   67.28 4.78 4.63   59.26   75.14 1.00     7549     7593
##    y0[11]   25.79   25.76 4.76 4.67   18.08   33.58 1.00     7490     7807
##    y0[12]   73.25   73.35 6.07 5.88   63.44   82.92 1.00     6285     6965
##    y0[13]   70.68   70.71 4.94 4.83   62.67   78.80 1.00     7997     7496
##    y0[14]   32.22   32.17 5.25 5.13   23.52   40.78 1.00     8116     7492
##    y0[15]   23.53   23.57 4.63 4.55   15.94   30.95 1.00     7439     7574
##    y0[16]   73.92   73.96 5.28 5.22   65.27   82.64 1.00     7296     6870
##    y0[17]   30.30   30.21 4.83 4.74   22.44   38.11 1.00     7869     7133
##    y0[18]   24.88   24.88 4.70 4.64   17.23   32.47 1.00     8104     7953
##    y0[19]   69.88   69.89 4.90 4.79   61.82   77.94 1.00     8002     6998
##    y0[20]   54.25   54.26 5.12 5.13   46.03   62.82 1.00     7878     7733
##    y0[21]   13.79   13.78 5.04 4.98    5.55   22.20 1.00     7359     7637
##    y0[22]   36.01   36.05 4.59 4.48   28.42   43.43 1.00     7823     7522
##    y0[23]   21.80   21.83 4.64 4.54   14.24   29.38 1.00     5683     7161
##    y0[24]   67.70   67.68 4.91 4.78   59.74   75.84 1.00     8215     7063
##    y0[25]   16.72   16.71 4.99 4.93    8.62   24.87 1.00     8032     7907
##    y0[26]   37.77   37.75 4.62 4.52   30.31   45.46 1.00     8221     7875
##    y0[27]   17.30   17.27 4.86 4.68    9.44   25.41 1.00     7736     7325
##    y0[28]   53.93   53.93 5.08 4.92   45.61   62.31 1.00     6234     5673
##    y0[29]   39.71   39.70 4.63 4.53   32.02   47.32 1.00     7320     7737
##    y0[30]    1.49    1.41 6.05 6.02   -8.30   11.48 1.00     3045     5252
##    y0[31]   28.83   28.84 4.60 4.42   21.37   36.43 1.00     7363     7283
##    y0[32]   36.17   36.13 4.59 4.53   28.69   43.70 1.00     7819     7187
##    y0[33]   48.75   48.85 5.07 4.99   40.34   57.01 1.00     8080     7857
##    y0[34]   49.71   49.71 4.89 4.85   41.67   57.59 1.00     6906     7733
##    y0[35]   19.27   19.36 4.87 4.72   11.24   27.07 1.00     6420     7923
##    y0[36]   39.72   39.72 4.60 4.54   32.15   47.31 1.00     7721     7491
##    y0[37]   39.70   39.72 4.63 4.54   32.09   47.19 1.00     7999     7803
##    y0[38]   26.91   26.91 4.95 4.86   19.00   35.14 1.00     8000     7650
##    y0[39]   12.86   12.88 4.89 4.73    4.81   20.93 1.00     6750     7059
##    y0[40]   25.45   25.48 4.70 4.62   17.66   33.11 1.00     7943     8014
##    y0[41]   37.45   37.43 4.63 4.59   29.94   44.99 1.00     8094     7680
##    y0[42]   28.44   28.40 4.85 4.75   20.55   36.34 1.00     8191     7636
##    y0[43]   71.07   71.06 4.95 4.96   63.12   79.23 1.00     8041     7913
##    y0[44]   19.04   19.00 4.76 4.66   11.23   26.92 1.00     8238     7365
##    y0[45]   19.80   19.83 4.95 4.92   11.69   27.96 1.00     7790     7353
##    y0[46]   62.15   62.12 4.81 4.70   54.27   70.02 1.00     7777     7664
##    y0[47]   91.75   91.74 5.60 5.58   82.48  100.86 1.00     7896     7053
##    y0[48]   33.02   33.04 4.57 4.47   25.54   40.55 1.00     7887     7521
##    y0[49]   63.82   63.89 4.81 4.67   55.84   71.64 1.00     7305     7229
##    y0[50]   40.47   40.44 4.70 4.60   32.77   48.19 1.00     8113     7939

m0=mcmc$draws('m0')
smm0=summary(m0)

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

xy=tibble(x=x,c=c,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)

qplot(x,y,shape=c,size=I(2),col=I('red'))+
  scale_shape_manual(values=1:k) + 
  #geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  #geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  #geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  #geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

qplot(x,y,facets=~c,size=I(2),col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

generalized linear mixed model

(X,y)i=1-n
b[b0,b1,...]
linear model    y~N(Xb,s)  
generalized linear model    y~dist.(m=link(Xb),s)  

fixed effect    b0, b1,...  
individual random effect   b0+r0i r0~N(0,sr0), b1+r1i r1~N(0,sr1),...  
class c=1-k  
class effect    b0+r0c r0~N(0,sr0), b1+r1c r1~N(0,sr1),...  

for y=dist.(m,s)
random intercept model    m=b0+r0i+b1*x, m=b0+r0c+b1*x  
random coefficient model  m=b0+(b1+r1i)*x, m=b0+(b1+r1c)*x  
mixed model   m=b0+r0i+(b1+r1i)*x, m=b0+r0c+(b1+r1c)*x  

note  
@ yi=b0+b1*xi+r0i is not useful, ri is included in s  
@ yi=b0+b1*xi+r0c is useful, repeated observation can be treated by class effect
@ when appling Poisson dist.(y is larger than 0, error is larger at large x),
  but variance is larger than mean (over dispersion),
  random intercept model is useful
n=20
x=runif(n,0,10)
r0=rnorm(n,0,1)
y=rpois(n,exp(2+0.1*x+r0))
qplot(x,y)

ex13-0.stan

generalized linear model, poisson regression

data{
  int N;
  vector[N] x;
  array[N] int y;
}
parameters{
  real b0;
  real b1;
}
model{
  y~poisson_log(b0+b1*x);
}
generated quantities{
  vector[N] m0;
  array[N] int y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=poisson_log_rng(m0[i]);
  }
}
mdl=cmdstan_model('./ex13-0.stan') 
data=list(N=n,x=x,y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = -4098.1 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       12        1479.7   0.000581323    0.00909063           1           1       19    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__    1479.70
##    b0         2.09
##    b1         0.21
##    m0[1]      2.31
##    m0[2]      2.89
##    m0[3]      3.45
##    m0[4]      2.70
##    m0[5]      3.41
##    m0[6]      2.71
##    m0[7]      3.67
##    m0[8]      3.78
##    m0[9]      2.89
##    m0[10]     3.64
##    m0[11]     3.77
##    m0[12]     3.85
##    m0[13]     3.78
##    m0[14]     3.77
##    m0[15]     3.13
##    m0[16]     3.48
##    m0[17]     2.67
##    m0[18]     3.08
##    m0[19]     2.46
##    m0[20]     3.97
##    y0[1]      5.00
##    y0[2]     13.00
##    y0[3]     55.00
##    y0[4]     22.00
##    y0[5]     34.00
##    y0[6]     22.00
##    y0[7]     48.00
##    y0[8]     42.00
##    y0[9]     17.00
##    y0[10]    29.00
##    y0[11]    45.00
##    y0[12]    52.00
##    y0[13]    43.00
##    y0[14]    42.00
##    y0[15]    18.00
##    y0[16]    35.00
##    y0[17]    19.00
##    y0[18]    17.00
##    y0[19]     4.00
##    y0[20]    59.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,'m')
##  variable    mean  median   sd  mad      q5     q95 rhat ess_bulk ess_tail
##    lp__   1478.72 1479.05 1.00 0.68 1476.68 1479.65 1.01      598      762
##    b0        2.10    2.10 0.14 0.13    1.88    2.35 1.02      228      209
##    b1        0.21    0.21 0.02 0.02    0.17    0.24 1.02      236      249
##    m0[1]     2.32    2.32 0.12 0.11    2.13    2.53 1.02      230      226
##    m0[2]     2.89    2.89 0.07 0.07    2.78    3.01 1.01      266      291
##    m0[3]     3.45    3.46 0.04 0.04    3.39    3.52 1.00     1271     1359
##    m0[4]     2.71    2.71 0.08 0.08    2.58    2.86 1.02      246      220
##    m0[5]     3.41    3.41 0.04 0.04    3.34    3.48 1.00      971     1379
##    m0[6]     2.72    2.72 0.08 0.08    2.58    2.86 1.02      246      220
##    m0[7]     3.67    3.67 0.04 0.04    3.60    3.75 1.00     1755     1523
##    m0[8]     3.78    3.78 0.05 0.05    3.69    3.86 1.00     1219     1598
##    m0[9]     2.90    2.90 0.07 0.07    2.79    3.02 1.01      267      291
##    m0[10]    3.64    3.64 0.04 0.04    3.57    3.71 1.00     1890     1533
##    m0[11]    3.77    3.77 0.05 0.05    3.69    3.85 1.00     1259     1598
##    m0[12]    3.85    3.85 0.05 0.05    3.76    3.93 1.00      866     1185
##    m0[13]    3.78    3.78 0.05 0.05    3.70    3.86 1.00     1176     1546
##    m0[14]    3.77    3.77 0.05 0.05    3.68    3.85 1.00     1275     1602
##    m0[15]    3.13    3.13 0.05 0.05    3.05    3.22 1.01      348      498
##    m0[16]    3.49    3.49 0.04 0.04    3.42    3.55 1.00     1447     1360
##    m0[17]    2.68    2.68 0.09 0.08    2.54    2.83 1.02      244      231
##    m0[18]    3.08    3.08 0.06 0.05    2.99    3.18 1.01      319      367
##    m0[19]    2.47    2.47 0.10 0.10    2.31    2.66 1.02      234      239
##    m0[20]    3.97    3.97 0.06 0.06    3.87    4.07 1.01      585     1000
##    y0[1]    10.28   10.00 3.52 2.97    5.00   16.00 1.00     1071     1822
##    y0[2]    18.13   18.00 4.47 4.45   11.00   26.00 1.00     1632     1856
##    y0[3]    31.50   31.00 5.83 5.93   22.00   41.00 1.00     2029     1947
##    y0[4]    15.26   15.00 4.10 4.45    9.00   22.00 1.00     1506     1754
##    y0[5]    30.32   30.00 5.57 5.93   22.00   40.00 1.00     1545     1886
##    y0[6]    15.09   15.00 4.15 4.45    9.00   22.00 1.00     1176     1574
##    y0[7]    39.56   39.00 6.76 7.41   29.00   51.00 1.00     2109     2105
##    y0[8]    43.67   43.00 7.09 7.41   32.00   56.00 1.00     1739     1715
##    y0[9]    18.21   18.00 4.34 4.45   11.00   26.00 1.00     1724     1839
##    y0[10]   38.11   38.00 6.09 5.93   29.00   48.00 1.00     2009     1629
##    y0[11]   43.51   43.00 6.93 7.41   33.00   55.00 1.00     2154     2016
##    y0[12]   46.92   47.00 7.44 7.41   35.00   60.00 1.00     1979     1889
##    y0[13]   43.90   44.00 6.98 7.41   33.00   56.00 1.00     1948     2023
##    y0[14]   43.25   43.00 7.11 7.41   32.00   55.00 1.00     1826     1931
##    y0[15]   22.94   23.00 4.81 4.45   15.00   31.00 1.00     1659     1822
##    y0[16]   32.55   32.00 5.96 5.93   23.00   43.00 1.00     2112     1810
##    y0[17]   14.53   14.00 3.96 4.45    8.00   22.00 1.00     1342     1960
##    y0[18]   21.74   22.00 4.75 4.45   14.00   30.00 1.00     1684     1927
##    y0[19]   12.00   12.00 3.70 2.97    6.00   18.00 1.00     1549     1894
##    y0[20]   52.72   53.00 7.91 7.41   40.00   66.00 1.00     1364     1839

m0=mcmc$draws('m0')
smm0=summary(m0)

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

xy=tibble(x=x,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)

qplot(x,y,size=I(2),col=I('red'))+
  scale_shape_manual(values=1:k) + 
  geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=exp(m)),xy,col='black')

ex13-1.stan

generalized linear mixed model, poisson regression

data{
  int N;
  vector[N] x;
  array[N] int y;
}
parameters{
  real b0;
  real b1;
  real<lower=0> sr0;
  vector[N] r0;
}
model{
  r0~normal(0,sr0);
  for(i in 1:N)
    y~poisson_log(b0+r0[i]+b1*x);
}
generated quantities{
  vector[N] m0;
  array[N] int y0;
  for(i in 1:N){
    m0[i]=b0+r0[i]+b1*x[i];
    y0[i]=poisson_log_rng(m0[i]);
  }
}
mdl=cmdstan_model('./ex13-1.stan') 
data=list(N=n,x=x,y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = 4594.93 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       99       29722.6      0.149664       4711.66      0.4869      0.7899      155    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      199       29874.8      0.193417   3.59653e+12      0.4627      0.4627      325    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      213       29877.3    0.00213885   2.78206e+11      0.6069      0.6069      343    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__   29877.30
##    b0         1.38
##    b1         0.32
##    sr0        0.00
##    r0[1]      0.00
##    r0[2]      0.00
##    r0[3]      0.00
##    r0[4]      0.00
##    r0[5]      0.00
##    r0[6]      0.00
##    r0[7]      0.00
##    r0[8]      0.00
##    r0[9]      0.00
##    r0[10]     0.00
##    r0[11]     0.00
##    r0[12]     0.00
##    r0[13]     0.00
##    r0[14]     0.00
##    r0[15]     0.00
##    r0[16]     0.00
##    r0[17]     0.00
##    r0[18]     0.00
##    r0[19]     0.00
##    r0[20]     0.00
##    m0[1]      1.71
##    m0[2]      2.61
##    m0[3]      3.49
##    m0[4]      2.33
##    m0[5]      3.42
##    m0[6]      2.33
##    m0[7]      3.83
##    m0[8]      3.99
##    m0[9]      2.62
##    m0[10]     3.78
##    m0[11]     3.99
##    m0[12]     4.11
##    m0[13]     4.00
##    m0[14]     3.98
##    m0[15]     2.99
##    m0[16]     3.54
##    m0[17]     2.28
##    m0[18]     2.90
##    m0[19]     1.96
##    m0[20]     4.30
##    y0[1]      7.00
##    y0[2]     14.00
##    y0[3]     35.00
##    y0[4]     13.00
##    y0[5]     32.00
##    y0[6]     13.00
##    y0[7]     50.00
##    y0[8]     56.00
##    y0[9]     18.00
##    y0[10]    40.00
##    y0[11]    65.00
##    y0[12]    51.00
##    y0[13]    46.00
##    y0[14]    45.00
##    y0[15]    11.00
##    y0[16]    37.00
##    y0[17]    13.00
##    y0[18]    11.00
##    y0[19]     5.00
##    y0[20]    84.00
mcmc=goMCMC(mdl,data,wrm=500)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 2 finished in 1.3 seconds.
## Chain 3 finished in 1.2 seconds.
## Chain 4 finished in 1.2 seconds.
## Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 1 finished in 1.5 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 1.3 seconds.
## Total execution time: 1.6 seconds.
seeMCMC(mcmc,'[m,r]')
##  variable     mean   median    sd   mad       q5      q95 rhat ess_bulk
##    lp__   29683.39 29682.80 15.22 16.01 29658.60 29709.90 1.11       38
##    b0         2.09     2.09  0.03  0.03     2.04     2.15 1.01      307
##    b1         0.21     0.21  0.00  0.00     0.20     0.22 1.01      324
##    sr0        0.01     0.01  0.01  0.00     0.00     0.02 1.11       34
##    r0[1]      0.00     0.00  0.01  0.00    -0.01     0.01 1.03     2958
##    r0[2]      0.00     0.00  0.01  0.00    -0.01     0.01 1.05     3644
##    r0[3]      0.00     0.00  0.01  0.00    -0.01     0.01 1.04     3359
##    r0[4]      0.00     0.00  0.01  0.00    -0.01     0.01 1.05     3047
##    r0[5]      0.00     0.00  0.01  0.00    -0.01     0.01 1.05     3155
##    r0[6]      0.00     0.00  0.01  0.00    -0.01     0.01 1.04     3246
##    r0[7]      0.00     0.00  0.01  0.00    -0.01     0.01 1.03     3608
##    r0[8]      0.00     0.00  0.01  0.00    -0.01     0.01 1.05     3414
##    r0[9]      0.00     0.00  0.01  0.00    -0.01     0.01 1.04     2472
##    r0[10]     0.00     0.00  0.01  0.00    -0.01     0.01 1.04     3250
##    r0[11]     0.00     0.00  0.01  0.00    -0.01     0.01 1.04     3013
##    r0[12]     0.00     0.00  0.01  0.00    -0.01     0.01 1.04     3278
##    r0[13]     0.00     0.00  0.01  0.00    -0.01     0.01 1.02     3034
##    r0[14]     0.00     0.00  0.01  0.00    -0.01     0.01 1.03     2981
##    r0[15]     0.00     0.00  0.01  0.00    -0.01     0.01 1.04     2324
##    r0[16]     0.00     0.00  0.01  0.00    -0.01     0.01 1.05     3144
##    r0[17]     0.00     0.00  0.01  0.00    -0.01     0.01 1.04     3057
##    r0[18]     0.00     0.00  0.01  0.00    -0.01     0.01 1.04     3203
##    r0[19]     0.00     0.00  0.01  0.00    -0.01     0.01 1.04     3609
##    r0[20]     0.00     0.00  0.01  0.00    -0.01     0.01 1.05     3338
##    m0[1]      2.31     2.31  0.03  0.03     2.26     2.36 1.01      344
##    m0[2]      2.89     2.89  0.02  0.02     2.86     2.92 1.01      408
##    m0[3]      3.45     3.45  0.01  0.01     3.43     3.47 1.00     1879
##    m0[4]      2.70     2.70  0.02  0.02     2.67     2.74 1.01      362
##    m0[5]      3.41     3.41  0.01  0.01     3.39     3.43 1.00     1555
##    m0[6]      2.71     2.71  0.02  0.02     2.67     2.75 1.00      365
##    m0[7]      3.67     3.67  0.01  0.01     3.65     3.69 1.00     2375
##    m0[8]      3.78     3.78  0.01  0.01     3.75     3.80 1.00     1994
##    m0[9]      2.89     2.89  0.02  0.02     2.86     2.92 1.01      407
##    m0[10]     3.64     3.64  0.01  0.01     3.62     3.66 1.00     2545
##    m0[11]     3.77     3.77  0.01  0.01     3.75     3.79 1.01     1670
##    m0[12]     3.85     3.85  0.01  0.01     3.83     3.87 1.00     1390
##    m0[13]     3.78     3.78  0.01  0.01     3.76     3.81 1.00     1787
##    m0[14]     3.77     3.77  0.01  0.01     3.75     3.79 1.00     1660
##    m0[15]     3.13     3.13  0.01  0.01     3.11     3.15 1.00      545
##    m0[16]     3.48     3.48  0.01  0.01     3.47     3.50 1.00     2215
##    m0[17]     2.67     2.67  0.02  0.02     2.64     2.71 1.01      372
##    m0[18]     3.08     3.08  0.02  0.02     3.05     3.10 1.00      542
##    m0[19]     2.47     2.47  0.03  0.03     2.42     2.51 1.00      332
##    m0[20]     3.97     3.97  0.02  0.02     3.95     4.00 1.00      928
##    y0[1]     10.09    10.00  3.23  2.97     5.00    16.00 1.00     1957
##    y0[2]     17.96    18.00  4.34  4.45    11.00    25.00 1.00     1712
##    y0[3]     31.46    31.00  5.79  5.93    22.00    41.05 1.00     2015
##    y0[4]     14.91    15.00  4.00  4.45     9.00    22.00 1.00     1896
##    y0[5]     30.26    30.00  5.54  5.93    21.00    39.05 1.00     1897
##    y0[6]     15.03    15.00  3.83  4.45     9.00    21.05 1.00     1874
##    y0[7]     39.36    39.00  6.12  5.93    30.00    50.00 1.00     1970
##    y0[8]     43.48    43.00  6.58  5.93    33.00    55.00 1.00     2050
##    y0[9]     17.92    18.00  4.19  4.45    11.00    25.00 1.00     1884
##    y0[10]    38.13    38.00  6.48  5.93    28.00    49.00 1.00     1990
##    y0[11]    43.38    43.00  6.68  7.41    32.00    55.00 1.00     1743
##    y0[12]    46.89    47.00  6.94  7.41    36.00    59.00 1.00     1918
##    y0[13]    43.97    44.00  6.65  5.93    33.00    55.00 1.00     1974
##    y0[14]    43.39    43.00  6.58  5.93    33.00    54.00 1.00     2138
##    y0[15]    23.09    23.00  4.85  4.45    15.00    31.00 1.00     2040
##    y0[16]    32.45    32.00  5.71  5.93    23.00    42.00 1.00     2013
##    y0[17]    14.56    14.00  3.80  4.45     8.00    21.00 1.00     1949
##    y0[18]    21.53    21.00  4.71  4.45    14.00    30.00 1.00     2139
##    y0[19]    11.71    11.00  3.38  2.97     7.00    17.05 1.00     1909
##    y0[20]    53.23    53.00  7.26  7.41    42.00    65.00 1.00     1886
##  ess_tail
##        78
##       377
##       329
##        56
##       581
##       598
##       687
##       765
##       600
##       906
##       527
##       478
##       517
##       605
##       735
##       554
##       527
##       771
##       580
##       824
##       557
##       663
##       656
##       678
##       474
##       686
##       755
##       587
##      1012
##       607
##      1018
##       793
##       993
##       980
##      1200
##      1139
##      1278
##      1386
##      1080
##      1039
##       583
##      1074
##       477
##       833
##      1941
##      1704
##      1675
##      1771
##      1936
##      1774
##      1953
##      1906
##      1895
##      1994
##      1799
##      1905
##      2025
##      2026
##      1985
##      1882
##      2050
##      1912
##      1858
##      1912

m0=mcmc$draws('m0')
smm0=summary(m0)

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

xy=tibble(x=x,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)

qplot(x,y,size=I(2),col=I('red'))+
  scale_shape_manual(values=1:k) + 
  geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=exp(m)),xy,col='black')

categorical dist. with parameters following Dirichlet dist.

h~Dir(a), h[0,1], sum(h)=1, a[0,1], sum(a)=1 
y~Cat(h), y[0,1], sum(y)=1

ex15-1.stan

data {
  int N;
  int K;
  array[N] int<lower=1,upper=K> y;
}
parameters {
  simplex[K] a;
  simplex[K] h;
}
model {
  h~dirichlet(a);
  y~categorical(h);
}
library(gtools)
n=20
h=rdirichlet(n,c(0.5,0.3,0.2))
summary(h)
##        V1              V2              V3       
##  Min.   :0.001   Min.   :0.000   Min.   :0.000  
##  1st Qu.:0.148   1st Qu.:0.032   1st Qu.:0.003  
##  Median :0.425   Median :0.082   Median :0.080  
##  Mean   :0.469   Mean   :0.295   Mean   :0.235  
##  3rd Qu.:0.791   3rd Qu.:0.696   3rd Qu.:0.313  
##  Max.   :0.996   Max.   :0.999   Max.   :0.927
c0=c(1,2,3)
for(i in 1:n) y[i]=sample(c0,1,prob=h[i,])
table(y) |> prop.table()
## y
##    1    2    3 
## 0.35 0.35 0.30
data=list(N=n,K=ncol(h),y=y)

mdl=cmdstan_model('./ex15-1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -26.0962 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       10      -22.6736   0.000289964   3.89817e-05      0.9996      0.9996       14    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -22.67
##      a[1]     0.34
##      a[2]     0.34
##      a[3]     0.32
##      h[1]     0.35
##      h[2]     0.35
##      h[3]     0.30
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -31.41 -31.07 1.48 1.31 -34.29 -29.66 1.00      915     1217
##      a[1]   0.34   0.32 0.18 0.19   0.08   0.66 1.00     1171      883
##      a[2]   0.34   0.32 0.18 0.19   0.08   0.66 1.00     1130      960
##      a[3]   0.32   0.30 0.17 0.18   0.08   0.63 1.00     1167     1244
##      h[1]   0.35   0.34 0.10 0.10   0.19   0.53 1.00     2774     1644
##      h[2]   0.35   0.35 0.10 0.11   0.19   0.53 1.00     2285     1639
##      h[3]   0.30   0.29 0.10 0.10   0.15   0.47 1.00     2261     1375

beta binomial dist.

y~betaB(n,a,b): y~B(n,p), p~beta(a,b)

ex15-2.stan

data {
  int N;
  int n;
  array[N] int y;
}
parameters {
  real<lower=0> a;
  real<lower=0> b;
}
model {
  a~cauchy(0,5);
  b~cauchy(0,5);
  y~beta_binomial(n,a,b);
}
generated quantities{
  int y1;
  y1=beta_binomial_rng(n,a,b);
}
n=30
a=3
b=4
p=rbeta(n,a,b)
n0=10
y=rbinom(n,n0,p)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    2.25    5.00    4.43    6.00    8.00
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,n=n0,y=y)

mdl=cmdstan_model('./ex15-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -246.925 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        8      -197.026    0.00201419    0.00162191           1           1       10    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__  -197.03
##      a        2.10
##      b        2.69
##      y1       3.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable    mean  median   sd  mad      q5     q95 rhat ess_bulk ess_tail
##      lp__ -196.14 -195.76 1.19 0.77 -198.33 -195.08 1.00      577      465
##      a       3.21    2.79 1.90 1.14    1.47    6.02 1.01      369      284
##      b       4.09    3.58 2.42 1.44    1.89    7.38 1.01      378      304
##      y1      4.39    4.00 2.38 2.97    1.00    8.00 1.00     1886     1943

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1')     

bimodal distribution, mixed normal distribution

n=100
y0=rnorm(n,0,1)
y1=rnorm(n,-5,1)
y2=rnorm(n,5,1)
y=sample(c(y0,y1,y2),n)
density(y) |> plot()

EM algorithm

library(mclust)

rst=Mclust(y)
summary(rst)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust E (univariate, equal variance) model with 3 components: 
## 
##  log-likelihood   n df  BIC  ICL
##            -239 100  6 -506 -507
## 
## Clustering table:
##  1  2  3 
## 32 35 33
rst$parameters
## $pro
## [1] 0.319 0.352 0.330
## 
## $mean
##       1       2       3 
## -4.7988 -0.0713  5.0721 
## 
## $variance
## $variance$modelName
## [1] "E"
## 
## $variance$d
## [1] 1
## 
## $variance$G
## [1] 3
## 
## $variance$sigmasq
## [1] 0.796
## 
## 
## $Vinv
## NULL
plot(rst)

ex17-1.stan

data {
  int N;
  vector[N] y;;
}
parameters {
  simplex[3] h; //ratio of mix
  real m1;
  real m2;
  real m3;
  real<lower=0> s1;
  real<lower=0> s2;
  real<lower=0> s3;
}
model {
  s1~cauchy(0,5);
  s2~cauchy(0,5);
  s3~cauchy(0,5);
  for (i in 1:N) {
    vector[3] p;
    p[1]=log(h[1]) + normal_lpdf(y[i] | m1, s1);
    p[2]=log(h[2]) + normal_lpdf(y[i] | m2, s2);
    p[3]=log(h[3]) + normal_lpdf(y[i] | m3, s3);
    target+=log_sum_exp(p);
  }
}
mdl=cmdstan_model('./ex17-1.stan')

data=list(N=n,y=y)
mle=mdl$optimize(data=data)
## Initial log joint probability = -375.528 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       99      -275.551   0.000193854     0.0354569           1           1      126    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      102      -275.551   0.000199822     0.0235199           1           1      129    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__  -275.55
##      h[1]     0.07
##      h[2]     0.89
##      h[3]     0.04
##      m1       0.35
##      m2       0.13
##      m3      -0.52
##      s1       0.17
##      s2       4.30
##      s3       0.02
mcmc=goMCMC(mdl,data,smp=2000)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:    1 / 2100 [  0%]  (Warmup) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:    1 / 2100 [  0%]  (Warmup) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:    1 / 2100 [  0%]  (Warmup) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:    1 / 2100 [  0%]  (Warmup) 
## Chain 4 Iteration:  101 / 2100 [  4%]  (Sampling) 
## Chain 3 Iteration:  101 / 2100 [  4%]  (Sampling) 
## Chain 4 Iteration: 1100 / 2100 [ 52%]  (Sampling) 
## Chain 3 Iteration: 1100 / 2100 [ 52%]  (Sampling) 
## Chain 2 Iteration:  101 / 2100 [  4%]  (Sampling) 
## Chain 4 Iteration: 2100 / 2100 [100%]  (Sampling) 
## Chain 4 finished in 1.4 seconds.
## Chain 1 Iteration:  101 / 2100 [  4%]  (Sampling) 
## Chain 3 Iteration: 2100 / 2100 [100%]  (Sampling) 
## Chain 3 finished in 1.6 seconds.
## Chain 2 Iteration: 1100 / 2100 [ 52%]  (Sampling) 
## Chain 2 Iteration: 2100 / 2100 [100%]  (Sampling) 
## Chain 2 finished in 2.1 seconds.
## Chain 1 Iteration: 1100 / 2100 [ 52%]  (Sampling) 
## Chain 1 Iteration: 2100 / 2100 [100%]  (Sampling) 
## Chain 1 finished in 2.5 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 1.9 seconds.
## Total execution time: 2.5 seconds.
seeMCMC(mcmc,ch=F)
##  variable    mean  median   sd  mad      q5     q95 rhat ess_bulk ess_tail
##      lp__ -246.99 -246.55 2.85 2.03 -251.08 -243.98 1.00     2453     1849
##      h[1]    0.34    0.33 0.05 0.05    0.26    0.42 1.03       85     1919
##      h[2]    0.33    0.33 0.05 0.05    0.25    0.41 1.03      101      579
##      h[3]    0.33    0.33 0.05 0.05    0.25    0.41 1.03      102     3823
##      m1      0.11   -0.05 3.99 2.61   -4.92    5.21 2.41        4       26
##      m2      1.34    3.55 4.10 2.86   -4.94    5.28 2.10        5       29
##      m3     -1.13   -2.04 4.07 4.06   -5.03    5.21 2.33        5       27
##      s1      0.99    0.91 1.48 0.16    0.71    1.27 1.16       16       89
##      s2      0.95    0.92 0.21 0.15    0.73    1.25 1.15       17       71
##      s3      0.99    0.96 0.19 0.17    0.73    1.32 1.22       12       76